# ===============================================
# Paris ratchet paper
# Code to run main models
# ===============================================

#### Table of contents ####

#' 1. Load data for models
#'  -- code to make these objects in "00_make_analysis_data.R"
#' 2. Run main models
#' 3. Full regression tables for appendix


#### Initialize #### 

library(tidyverse)
library(broom)
library(modelsummary)


## Check working directory
list.files() 

#' Should see:
#' [1] "_readme_ratchet.html" "_readme_ratchet.Rmd"  "code"
#' [4] "grateful-refs.bib"    "input"                "output"
#' [7] "plots"  "text"
#' 
#' If not, set working directory properly -- try
# setwd("..")

# Define core ISO3C codes to use in analysis
core_iso3c <- c("AFG", "AGO", "ALB", "AND", "ARE", "ARG", "ARM", "ATG", "AUS", "AUT", "AZE", "BDI", "BEL", "BEN", "BFA", "BGD", "BGR", "BHR", "BHS", "BIH", "BLR", "BLZ", "BOL", "BRA", "BRB", "BRN", "BTN", "BWA", "CAF", "CAN", "CHE", "CHL", "CHN", "CIV", "CMR", "COD", "COG", "COK", "COL", "COM", "CPV", "CRI", "CUB", "CYP", "CZE", "DEU", "DJI", "DMA", "DNK", "DOM", "DZA", "ECU", "EGY", "ERI", "ESP", "EST", "ETH", "FIN", "FJI", "FRA", "FSM", "GAB", "GBR", "GEO", "GHA", "GIN", "GMB", "GNB", "GNQ", "GRC", "GRD", "GTM", "GUY", "HND", "HRV", "HTI", "HUN", "IDN", "IND", "IRL", "IRN", "IRQ", "ISL", "ISR", "ITA", "JAM", "JOR", "JPN", "KAZ", "KEN", "KGZ", "KHM", "KIR", "KNA", "KOR", "KWT", "LAO", "LBN", "LBR", "LBY", "LCA", "LIE", "LKA", "LSO", "LTU", "LUX", "LVA", "MAR", "MDA", "MDG", "MDV", "MEX", "MHL", "MKD", "MLI", "MLT", "MMR", "MNE", "MNG", "MOZ", "MRT", "MUS", "MWI", "MYS", "NAM", "NER", "NGA", "NIC", "NLD", "NOR", "NPL", "NRU", "NZL", "OMN", "PAK", "PAN", "PER", "PHL", "PLW", "PNG", "POL", "PRK", "PRT", "PRY", "QAT", "ROU", "RUS", "RWA", "SAU", "SDN", "SEN", "SGP", "SLB", "SLE", "SLV", "SOM", "SRB", "SSD", "STP", "SUR", "SVK", "SVN", "SWE", "SWZ", "SYC", "SYR", "TCD", "TGO", "THA", "TJK", "TKM", "TLS", "TON", "TTO", "TUN", "TUR", "TUV", "TZA", "UGA", "UKR", "URY", "USA", "UZB", "VCT", "VEN", "VNM", "VUT", "WSM", "YEM", "ZAF", "ZMB", "ZWE")

### Load data

data_for_models <- readRDS("output/data_for_models20230920.Rds")

## Table settings

coefficients_renamed = c(
  "io_percentage" = "IO Paris",
  "tradeflow_percentage" = "Trade Paris",
  "eite_percentage" = "Trade (EITE) Paris",
  "green_percentage" = "Trade (Green) Paris",
  "competition_tradeflow_percentage" = "Trade (competition) Paris",
  "competition_eite_percentage" = "Trade (competition, EITE) Paris",
  "trade_openness" = "Trade openness",
  "renewables_share" = "Renewable electricity",
  "industry" = "Industry",
  "fossil_rents" = "Fossil rents",
  "paris_target_safe" = "Paris target",
  "(Intercept)" = "(Intercept)")

gm <- tibble::tribble(
  ~raw, ~clean, ~fmt,
  "nobs", "Observations", 0,
  "r.squared", "R2", 2)


# Function to detect whether controls are in a model
detect_controls <- function(model, control_vars) {
  model_terms <- names(coef(model))
  has_controls <- any(control_vars %in% model_terms)
  return(ifelse(has_controls, "Yes", "No"))
}

# Define a control variable to look for
control_vars_flag <- c("fossil_rents")


#### Table 1, trade ties ####

t1_a <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + tradeflow_percentage,
     data = .)
t1_b <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + tradeflow_percentage +
       trade_openness + industry + renewables_share + fossil_rents,
     data = .)
t1_c <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + eite_percentage,
     data = .)
t1_d <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + eite_percentage +
       trade_openness + industry + renewables_share + fossil_rents,
     data = .)
t1_e <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + competition_tradeflow_percentage,
     data = .)
t1_f <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + competition_tradeflow_percentage +
       trade_openness + industry + renewables_share + fossil_rents, 
     data = .)
t1_g <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + competition_eite_percentage,
     data = .)
t1_h <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + competition_eite_percentage +
       trade_openness + industry + renewables_share + fossil_rents,
     data = .)

## Look for controls in each model 
t1_controls_status <- c(
  t1_a = detect_controls(t1_a, control_vars_flag),
  t1_b = detect_controls(t1_b, control_vars_flag),
  t1_c = detect_controls(t1_c, control_vars_flag),
  t1_d = detect_controls(t1_d, control_vars_flag),
  t1_e = detect_controls(t1_e, control_vars_flag),
  t1_f = detect_controls(t1_f, control_vars_flag),
  t1_g = detect_controls(t1_g, control_vars_flag),
  t1_h = detect_controls(t1_h, control_vars_flag)
)

## Convert to tibble and add attributes to slot into modelsummary()
t1_latex_controls <- tibble::tibble(term = "Controls", !!!t1_controls_status)
attr(t1_latex_controls, 'position') <- c(13)


table1_latex <- modelsummary(models = list(t1_a, t1_b, t1_c, t1_d, 
                           t1_e, t1_f, t1_g, t1_h),
             output = "latex_tabular",
             fmt = 2,
             coef_map = coefficients_renamed,
             coef_omit = c("trade_openness|industry|renewables_share|fossil_rents"),
             gof_map = gm, 
             add_rows = t1_latex_controls,
             # notes = list("Outcome is Glasgow mitigation target"),
             stars = c("+" = 0.1, "*" = 0.05, "**" = 0.01))
table1_latex
#' -> paste into manuscript
#' Add "\midrule" before controls 
#' Bold the spatial matrix term
#' Math the R2 term
#' Check caption



#### Table 2, Political ties ####

t2_a <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + io_percentage,
     data = .)
t2_b <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + io_percentage +
       trade_openness + industry + renewables_share + fossil_rents,
     data = .)
t2_c <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + io_percentage + tradeflow_percentage,
     data = .)
t2_d <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + io_percentage + tradeflow_percentage +
       trade_openness + industry + renewables_share + fossil_rents,
     data = .)
t2_e <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + io_percentage + eite_percentage,
     data = .)
t2_f <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + io_percentage + competition_tradeflow_percentage,
     data = .)
t2_g <- data_for_models %>%
  lm(glasgow_target_safe ~ paris_target_safe + io_percentage + competition_eite_percentage,
     data = .)


## Look for controls in each model 
t2_controls_status <- c(
  t2_a = detect_controls(t2_a, control_vars_flag),
  t2_b = detect_controls(t2_b, control_vars_flag),
  t2_c = detect_controls(t2_c, control_vars_flag),
  t2_d = detect_controls(t2_d, control_vars_flag),
  t2_e = detect_controls(t2_e, control_vars_flag),
  t2_f = detect_controls(t2_f, control_vars_flag),
  t2_g = detect_controls(t2_g, control_vars_flag)
)

## Convert to tibble and add attributes to slot into modelsummary()
t2_latex_controls <- tibble::tibble(term = "Controls", !!!t2_controls_status)
attr(t2_latex_controls, 'position') <- c(15)

## Make table
table2_latex <- modelsummary(models = list(t2_a, t2_b, t2_c, t2_d, 
                           t2_e, t2_f, t2_g),
             # output = "latex_tabular",
             fmt = 2,
             coef_map = coefficients_renamed,
             coef_omit = c("trade_openness|industry|renewables_share|fossil_rents"),
             gof_map = gm, 
             add_rows = t2_latex_controls,
             # notes = list("Outcome is Glasgow mitigation target"),
             stars = c("+" = 0.1, "*" = 0.05, "**" = 0.01))
table2_latex 
#' -> paste into manuscript
#' Add "midrule" above controls
#' Bold the spatial matrix term
#' Math the R2 term
#' Check caption

#### Full regression tables ####

## Convert to tibble and add attributes to slot into modelsummary()
t1_full_controls <- tibble::tibble(term = "Controls", !!!t1_controls_status)
attr(t1_full_controls, 'position') <- c(21)

## Make table without omitting covariates
table1_full <- modelsummary(models = list(t1_a, t1_b, t1_c, t1_d, 
                                         t1_e, t1_f, t1_g, t1_h),
                           output = "latex_tabular",
                           fmt = 2,
                           coef_map = coefficients_renamed,
                           # coef_omit = c("trade_openness|industry|renewables_share|fossil_rents"),
                           gof_map = gm, 
                           add_rows = t1_full_controls,
                           # notes = list("Outcome is Glasgow mitigation target"),
                           stars = c("+" = 0.1, "*" = 0.05, "**" = 0.01))
table1_full
#' -> paste into manuscript
#' Add "midrule" above controls
#' Bold the spatial matrix term
#' Math the R2 term
#' Check caption


## Convert to tibble and add attributes to slot into modelsummary()
t2_full_controls <- tibble::tibble(term = "Controls", !!!t2_controls_status)
attr(t2_full_controls, 'position') <- c(23)

table2_full <- modelsummary(models = list(t2_a, t2_b, t2_c, t2_d,
                                           t2_e, t2_f, t2_g),
                             output = "latex_tabular",
                             fmt = 2,
                             coef_map = coefficients_renamed,
                             # coef_omit = c("trade_openness|industry|renewables_share|fossil_rents"),
                             gof_map = gm, 
                             add_rows = t2_full_controls,
                             # notes = list("Outcome is Glasgow mitigation target"),
                             stars = c("+" = 0.1, "*" = 0.05, "**" = 0.01))
table2_full
#' -> paste into manuscript
#' Add "midrule" above controls
#' Bold the spatial matrix term
#' Math the R2 term
#' Check caption






